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Requirements of small-size and lightweight heat exchanger devices in power, space, 
process and aerospace industries have resulted in the development of specially designed 
heat transfer surfaces. However, in most practical cases, the heat transfer coefficients 
are low. In order to enhance heat transfer between the flowing fluid and closely spaced 
parallel plate channels, in the case of plate-fin heat exchangers, protrusions can be 
mounted on the channel surfaces. The protrusions in the form of winglet-type vortex 
generators are capable of etihancing heat transfer at the expense of relatively small 
increase in pressure penalty as compared to other protrusions geometries. 

To maximize heat transfer while incumng minimum pressure penalty is an important 
design goal and this is the primary objective of the present investigation. To evaluate 
different heat exchanger geometries this study utilizes the MAC (Marker and Cell) 
algorithm. A three-dimensional numerical model of this algorithm is used to solve Ml 
Navier-Stokes equations together with the governing equations for energy conservation 
in a rectangular channel with a built-m winglet vortex generator. 

For the bi-criteria optimization of heat exchanger design geometry, multiple criteria 
optimization methods may be used. This study has appraised the usefulness of EVOP 
for solving multi-objective engineering design problems for which GA (Genetic 
Algorithms) or other multi-objective methods cannot work due to practical limitations. 
EVOP is a simple but powerful technique, created by Box and Draper (1969), which has 
been widely and successMly used in empirical process optimization in a variety of 
process industries. 

This work has demonstrated that EVOP can be an effective method to find Pareto 
Optimal engineering designs particularly when the number of design evaluations 
performed must remain small in number. 


List of Figures 


Figure No. Title Page No. 

1 . 1 Productivity improvement with EV OP technique 3 

1 .2 Size of improvement and Noise level 4 

1 .3 Controlled Vanation in Design Inputs 5 

1 .4 The Efficient Front in Bi-objective Maximization Problem 1 0 

2. 1 Typical arrangement of heat exchanger cores 

(a) Gas-hquid fin-tube cross-flow 13 

(b) Plate fin (single or multi-pass) 13 

2.2 Various methods of making enhanced tubes: 

(a) helical rib on inner surface and integral fin on outer surface 1 5 

(b) internal fins on inner surface and porous coating on other surface 15 

(c) twisted tape insert on inner surface and integral fins on outer 15 
surface 

2.3 (a)Plate -fin heat exchanger and its surface geometries 16 

(b) plain rectangular fins 16 

(c) plain triangular fins - 16 

(d) wavy fins 16 

(e) offset strip , 16 

(f) perforated fins 16 

(g) louvered fins; after Webb (1987) 16 

2.4 Boundary conditions and fictitious boundary cells 20 

2.5 (a) Grid spacing in the computational domain and the location 22 

of the wing type vortex grid showing the descretized variables 
(b) Three dimensional staggered grid showing the locations 22 

of the descretized variables. 

2.6 Relative location of velocity components and wing on X-Y plane 24 

atZ=0 

2.7 Side edge of the wing-type vortex generator on Y-Z plane 25 



Figure No. 


Title 


Page No. 

2.8 Thermal boundary conditions on the wing-type vortex 25 

2.9 Periodic boundary conditions for pressure and velocities on the 26 

top and the bottom wall. 

3.1 Geometry of Channel 39 

3.2 Geometry of winglet 39 

3.3 (a) Symbolic Spanwise Average Nusselt Number Distribution 41 

(b) Heat Exchanger with Upper and Bottom Plates showing cells in 41 
flow direction. 

3.4 Symbolic representation of Bi-criteria problem solutions in the 42 

Response Space 

3.5 Two dimensional variants in the Design Space 43 

3.6 New Pareto front in the Response Space after new input runs 43 

4.1 Channel Geometry 

(a) Geometry of heat exchanger plates 46 

(b) Geometry of winglet 46 

(c) Top view of the winglet on heat-exchanger plate 46 

4.2 (a) Channel dimensions 47 

(b) Sample Winglet Dimensions 48 

4.3 EVOP Table for Initial Runs in the input space 50 

4.4 Initial Performance 5 1 

4.5 (a) EVOP table for second phase 52 

(b) Result after second phase of EV OP 53 

4.6 (a) Table for 3"*^ phase of EVOP 53 

(b) Results after 3’^*^ phase of EVOP 54 

4.7 Pareto Solutions 56 

4.8 EVOP Table for determination of Input sets of 4* phase 57 

4.9 Output of all the runs 58 

4.10 Pareto Optimaland Non-Pareto Design outputs 59 

4.1 1 EVOP Table for Pareto Optimal and Non-Pareto designs 60 



List of Tables 


Table No. 

4.1 

4.2 

4.3 

4.4 

4.5 
5.1 


Title 

Summary of results after third phase of 
EVOP 

List of Pareto Solutions after third phase 
Set of new inputs for fourth phase of EVOP 
Inputs and outputs of all Designs 
Inputs and Outputs of Pareto Optimal Designs 
Final Set of Pareto Optimal Designs 


Page No. 

55 

56 
58 
61 
62 
63 



1. Introduction 


1.1 A Heat Exchanger Design Problem 

Requirements of small-size and lightweight heat exchanger devices in power, process and 
aerospace industries have resulted in the development of specially designed heat transfer 
surfaces. It is possible to distinguish different basic flow configurations, such as internal 
flow through tubes, external flow normal to the tubes and channel flows in a wide range 
of heat exchangers. In most practical cases, the heat transfer coefficients for the 
configuration of flow over flat surfaces are significantly low. In order to enhance heat 
transfer between the flowing fluid and closely spaced parallel plate channels, in the case 
of plate-fin heat exchangers, protrusions can be moimted on the channel surfaces. The 
protrusions in the form of slender delta wing or winglet-type vortex generators are 
capable of enhancing heat transfer at the expense of relatively small increase in pressure 
penalty as compared to other protrusion geometries. 

Experimental investigations reveal enhancement of heat transfer in the presence of 
longitudinal vortices. Such vortices are characterized by their rotating motion and 
longevity. (Biswas et al. 1995) 

Computational studies on related topics have been performed by Fiebig and Biswas 
(1989) and Chattopadhyay (1992) in a geometrical configuration of delta wing placed 
inside a rectangular channel. Both studies confirm significant transport enhancement. In 
another recent study, Biswas et a/. (1994) have determined that winglets are a more 
attractive choice than delta wings as vortex generators for enhancing heat transfer. 

The present study is aimed at a numerical investigation of heat transfer enhancement 
while incurring minimized pressure penalty in a channel with a built-in winglet type 
vortex generator placed in a fully developed laminar flow. Since this study utilizes the 
state-of-art MAC (Marker and Cell ) algorithm to do design evaluation with each 
evaluation taking 3-4 hours on a SGI Origin 200 (4 Mips RIOOOO processor, 512 MB 
Ram, 30 GB HDD, 10/100 Ethernet Interface machine), GA- based multi-objective 
methods such as ENGA (Bagchi, 1999) could not be used here. Rather, we adopted a 
variation of the EVOP method described below for this purpose. 
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1.2 EVOP Technique for Optimization 


1.2.1 Theoretical Foundation of EVOP 

Over the course of its life a typical industrial process passes through many stages of 
development. When a plant is built, it ideally incorporates the best design possible, given 
the available knowledge and resources. Nevertheless the process of ‘tuning’ still remains 
to be done to improve quality, yield, cost, performance, cycle time etc. 

Consequently in most industrial organizations considerable effort is made to further 
improve the process after its start-up. In addition to the recording and examination of 
routine plant data, special studies are frequently carried out in the laboratory, on the pilot 
plant, and also on the full-scale process by the research and development staff. Such 
investigations, particularly when they employ efficient modem investigatory tools of 
statistical design and analysis, are most valuable and must continue (Box and Draper, 
1969). 

The shortage of technical personnel however inevitably limits the amount of special 
investigation of this kind. There is a special technique that makes only sparing use of 
technical manpower and complements the special investigations already referred to. 
Fig.(l.l) shows the typical improvement in productivity with the use of EVOP technique 
as compared to the productivity improvement with normal technical effort alone. (Box 
and Draper, 1969). 

Because of the continual technical effort exerted on process improvement, we normally 
find a steady increase in productivity. From one point of view, one does not know at the 
beginning of any particular time which particular processes will be improved or what the 
precise nature of the improvements will be but, nevertheless, one can predict with some 
exactitude what the overall effect of applying his technical work will be. From another 
point of view, the near certainty of the specific amount of process improvement expected 
during the forthcoming year may be depressing. One knows that without some new 
initiative he cannot hope to increase the rate of improvement of productivity. 

There is a simple but powerful empirical technique, that extends one such initiative. 
Created by Box and Draper (1969) it is called Evolutionary Operation or more briefly. 
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EVOP. It has now been widely and successfully used in the chemical industry, but would 
undoubtedly be of value in other types of manufacture as well. The present study is an 
attempt to use EVOP to find Pareto optimal heat exchanger designs. 

Evolutionary Operation is really a management technique in which a continuous routine 
becomes the basic mode of operation for the plant and replaces normal static operation. 
Evolutionary Operation is not a substitute for fundamental R & D investigation. Rather, 
such investigation should use full factorial experiments (DOE) and should continue side 
by side with EVOP. Evolutionary Operation does, however, often indicate areas to which 
more fundamental approaches such full-fledged DOE could usefully be directed. 



Year 

Figure 1.1: Productivity improvement with EVOP technique 


In the EVOP method a carefully planned cycle of minor variants on the process is agreed 
upon by those seeking to improve the process. The routine of EVOP consists of running 
each of the variants in turn and continually repeating the cycle. The cycle of variants 
follows a simple pattern, the persistent repetition of which allows evidence concerning 
yield and physical properties of the product in the immediate vicinity of the works 
process to accumulate during routine process runs. In this way we use routine process 
running to generate not only the product we require for business but also the information 
we need to improve it. 
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Controlled variations thus being introduced into the process, the effect of selection is 
introduced by arranging for the results to be continuously presented. The stream of 
information concerning the outputs from the various process conditions is summarized. 
The information is set out in such a way that the analyzer can at any time see what weight 
of evidence exists for moving the center of scheme of variants to some new point, what 
types of changes are undesirable from the standpoint of producing results of inferior 
quality, how much the scheme is costing to run, and so on. 

In making permanent change in the routine of process, the situation is very different from 
that which we meet in running specialized experiments such as those using DOE (Design 
of Experiments). The latter will last a limited time, during which special facilities can be 
made available. Evolutionary Operation, however, is virtually a permanent method of 
running the process. Changes in the levels of the variables can be permitted whose effects 
are virtually undetectable in individual mns, and only techniques simple enough to be run 
continuously by works personnel themselves under actual conditions of manufacture can 
be employed. 

Suppose a new plant has just been built. We imagine the improvements, which could 
arise from suitable adjustments to the new process measured on some scale such as 
profitability and arranged in a descending order of magnitude. (Refer Fig. 1 .2) 



Figure 1.2: Size of improvement and Noise level 
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A few very large effects to the left of the diagram will be quickly detected after startup. 
Their causes will then be tracked down and adjustments and improvement of the process 
will result. 

The reason that these effects can be rapidly detected, hence, exploited, is that the ‘signal’ 
they produce is large compared with the underlying level of variation of the process, 
sometimes called ‘noise level’. 

Noise arises from a variety of sources such as variabihty of raw materials, inability to 
maintain precisely input variables at set levels, and from instrument and measurement 
error. Such variation is often regarded as inherent. 

To discover effects buried in noise we must improve the signal-to-noise ratio. We must 
either decrease the effective noise level or increase the signal level. In EVOP we do both. 
The signal is deliberately increased by introducing changes of a carefully chosen kind in 
the process conditions or variables under study. Repetition of the changes and averaging 
of the results reduce the effective noise level. 

1.2.2 An Example of EVOP 



Concentration, % 

Figure 1.3: Controlled Variation in Design Inputs 
The variables under study 

At a particular stage of a development of a certain batch process, two process variables or 
factors were being studied. These were, respectively, the percentage concentration of a 
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certain feed material and the temperature at which the reaction was conducted. 


A Pattern of Variants of Design 

The scheme of variants is shown in above figure. The current works process is labeled 
‘0’ and the four variants are labeled 1,2,3, and 4. One batch of product was made at each 
set of conditions, which were run successively in the order 0,1, 2,3,4; 0,1, 2,3,4; and so on. 

Responses of Interest 

Three responses were recorded; 

1 . The cost of manufacturing unit weight of product. 

2. The percentage of a certain of impurity. 

3 . A measure of fluidity. 

The Information Board 

The information was displayed on a “board”. The essential thing is that this information 
(data) and its display should be simple. In general, the information should consist of cycle 
number, phase number and the output (response) factors in figures. 

Analysis of the Information Board 

After studying the information board and combining the ideas it conveys with production 
requirements and his special knowledge of the process, the superintendent can make one 
of two basic decisions: 

1. To allow the scheme to proceed unchanged and wait for additional information from 
the next cycle. - 

2. To modify operation in some way, so beginning another phase of the EVOP scheme. 
Under alternative 2, a number of possibilities are available. Some of these are as follows: 

• Adopt one of the evolved 

• variants as the new ‘works process’ and recommence the cycle about this new 
center point. 
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• Explore an indicated favorable direction of advance and recommence the 
cycle about the best condition found. 

• Substitute new variables for one or more of the old variables. 

• Modify the pattern of variables for one or more of the old variables waiting to 
be tested on the EVOP scheme. 

We can see from this example that the actions to be taken in various circumstances are 
not precisely specified. H\nnan judgement is a very important part of an EVOP 
investigation as it should be any investigation. A properly run EVOP scheme ensures that 
the process out of the hands of the process superintendent is continuously fed with clear 
information about this process on which he then takes what actions he deems appropriate. 
Although there are certain principles in this game, judgement is called for at every stage 
to evaluate the possibilities of future success or failure and to act accordingly. 

1.3 Multi-Objective Optimization 

Many real world problems involve multiple measures of performance, or objectives, 
which should be optimized simultaneously. In certain cases, objective functions may be 
optimized separately from each other and insight gained concerning the best that can be 
achieved in each performance dimension. However, suitable solutions to the overall 
problem can seldom be found in this way. Optimal performance according to one 
objective, if such an optimum exists, often implies unacceptably low performance in one 
or more of the other objective dimensions, creating the need for a compromise to be 
reached (Fonesca and Fleming, 1995). A suitable solution to such problems involving 
conflicting objectives should offer ‘acceptable’, though possibly sub-optimal in the 
single-objective sense, performance in all objective dimensions, where ‘acceptable’ is a 
problem-dependent and ultimately subjective concept. 

1.4 Multiple Criteria Decision Making 

Decision making is the selection of an act or courses of action from among alternative 
acts or courses of actions such that it will produce optimal results under some criteria of 
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optimization (Tabucanon, 1989). This concise definition of decision making invokes 
further elaboration to a certain extent. Before the problems can be considered well 
defined, the set of alternatives and the set of criteria have to be known and established 
first; only then can the selection process commence. What makes multiple criteria 
decision-making complex is the plurality of the criteria involved in the problem. In 
decision analysis of complex systems, such terms as ‘multiple criteria’, ‘multiple 
objectives’, or ‘multiple attributes’ are used interchangeably to describe decision 
situations while there are no universal definitions of these terms. Multiple criteria 
decision-making (MCDM) has seemed to emerge as the accepted nomenclature for all 
models and techniques dealing with multiple objective decision making or multiple 
attribute decision making. 

An example of the Multiple Criteria Decision Making 
Production Planning involves the following multiple objectives: 

Max {Total net revenue} 

Max {Minimum net revenue in any period} 

Min {Backorders} 

Min{Overtime} 

Min {Finished goods inventory} 

Such multiplicity of criteria occur in almost every area of business decision making and 
operations. They appear in accounting, finance, operations management, manpower 
planning etc. 

The general multi-objective problem requiring the optimization of k objectives may be 
formulated as follows: 


Max (min) Zj = ^(X), j=l,2,. . .,k 
Subject to 

• g,(X) <b„ i= l,2,...,m 

• X > 0 
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where X is a vector of decision variables and gi(X) are the inequality constraints. In 
multi-criteria optimization, because the objectives are often conflicting, there does not 
exist a single unique solution, which is globally optimum with respect to all the 
conflicting objectives. The increment in any of these objectives will decrease the others 
and vice versa. 

A necessary condition of MCDM is the presence of more than one criteria. The sufficient 
condition is that the cnteria must be conflicting in natirre. In summary, a problem can be 
considered as that of MCDM if and only if there appears at least two conflicting criteria 
and there are at least two alternative solutions (Tabucanon, 1989). 

Criteria are said to be in conflict if the full satisfaction of one will result in impairing or 
precluding the full satisfaction of the other(s). The criteria are considered to be ‘strictly’ 
conflicting if the increase in satisfaction of one results in a decrease in satisfaction of the 
other. The sufficient condition of MCDM, however, does not necessarily stipulate 
‘strictly’ conflicting criteria. 

In view of the conflicting nature of the criteria involved in MCDM, ciioosing the ‘best’ 
alternative is indeed a difficult task for the decision-maker. Consequently there is a need 
for methods to systematically resolve the conflicts among criterion (or objectives) in 
order to reach acceptable compromises and make up with satisfying (or termed as 
‘satisficing’ ) solutions. 


1.5 The Concept of Pareto Optimality 

A special “family” of solutions of a multi-objective optimization comprises all those 
elements of the search space- which are such that the components of the corresponding 
objectives vectors caimot be all simultaneously improved. In economics such solutions 
are called Pareto optimal. 

A more formal definition of Pareto optimality is as follows (Tabucanon, 1989): 

Consider, without loss of generality, the minimization of the n components fk,k= 1,. . .,n, 
of a vector variable x in a universe Ux is said to be Pareto optimal if and only if there is 
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no Xv e Ux for which v =fixv) = (v;_.. y„) dominates u = (ui, Un), i.e., there is no 

Xv e Ux such that 

Vie v,<M,n 3i€ {\,...,n}\vi<Ui 

The set of all Pareto-optimal decision vectors is called the Pareto optimal, efficient, or 
admissible set of the solutions. The corresponding set of objective vectors is called the 
non-dominated set. In practice, however, it is not xmusual for these terms to be used 
interchangeably to describe solutions of a multi-objective optimization problem. 

In Fig. 1.4, the Pareto optimal set of solutions is shown by dark points on the dotted line, 
which are solutions numbered 10, 1 1,8 and 7. 
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Figure 1.4;The Efficient Front in Bi-objective 
Maximization Problem 


The notion of Pareto Optimality is only a first step towards the practical solution of a 
multi-objective problem. This decision involves selecting subsequently a single 
compromise solution from the non-dominated set according to some preference 
information. All other solutions are dominated in the Pareto sense and they may be 
ignored. 
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Finding Pareto optimal solutions, however, is not straight forward. Recently, methods 
based on metaheuristics have been suggested (Srinivas and Deb, 1995; Bagchi, 1999). 
The present work is an attempt to use EVOP when such methods would not be 
practicable. 
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2. Problem Definition 


2.1 Description of the Problem 

2.1.1 Introduction to Heat Exchanger Design 

Fin-tubes are generally used in gas-liquid heat exchangers where the gas is in the 
crossflow to the tubes and the liquid flows inside the tubes. Figure (2.1) shows the 
arrangement of the core region of a fin-tube crossflow heat exchanger. The purpose of the 
fin is to enhance the heat transfer coefficient on the gas side to a value comparable to the 
liquid side. The area ratio of the fins to the tubes varies depending on the application and 
may reach a value of 30. Plate-fin heat exchangers also find application in many 
industrial processes. Such heat exchangers may be used to exchange energy between two 
different gases or liquids. Figure (2.1) illustrates the basic arrangement of the core region 
of a plate-fin heat exchanger. In forced convection heat transfer between a gas and a 
liquid, the heat transfer coefficient of the gas may be 10 to 50 times smaller than that of 
the liquid. In order to increase tihe heat transfer on the gas side, the fin geometry can be 
manipulated (Webb, 1987) in various ways. In addition, compactness of the heat 
exchangers is demanded for effective utiHzation of energy and to bring about a save on 
material and space. The performance of heat exchanger surfaces brings about a save on 
material and space. The performance of heat exchanger surfaces can be enhanced by 
mounting protrusions on the surfaces. The offset strip-fin and the louvered fins have been 
widely used for this purpose (Webb, 1987). The basic mechanism involved in the above 
mentioned cases is the periodic interruption of the boundary layer by separation, wake 
recovery and generation of thin developing boundary layer with high Nusselt Number. 

A somewhat different concept for the reduction of the thermal resistance is to induce 
longitudinal streamwise vortices in the flow field. The longitudinal vortices can be 
generated by mounting delta wing or winglet type vortex generators on the flat surfaces 
(Fig. 2.1). The longitudinal vortices develop along the side edge of the wing-shaped 
vortex generator due to file pressure difference between the front surface facing flow and 
the back surface. These streamwise vortices interact with an otherwise two-dimensional 
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boundary layer and produce a three dimensional swirling flow that mixes near wall fluid 
with the free stream. This mechanism strongly enhances the exchange of fluid between 



Figure (2.1): Typical arrangement of heat exchanger cores (a)Gas-liquid fin-tube cross- 
flow (b) Plate fin (single or multi pass) 

the wall and the core region of the flow field, which causes high heat transfer 
augmentation. The additional pressure lasses are modest because the form drag for such 
wing-type slender bodies is low. A complete numerical evaluation of heat transfer m a 
fin-tube or plate-fin crossflow heat exchanger would be extremely difficult. In order to 
understand the basic mechanisms involved in enhancement of heat transfer and to 
evaluate the performance parameters quantitatively, a need is felt to analyze three- 
dimensional flow and heat transfer in a horizontal channel with built-in vortex generators. 
A number of experimental investigations (Edwards and Alker, 1974; Russel et al, 1982, 
Fiebig et al, 1986) in fields relevant to the present problem have been reported in the 
literature. Computational studies on related topics have been performed by Fiebig et al 
(1989) and Biswas and Chattopadhyay (1992) for laminar flows. 
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2.1.2Enha]icement of Heat Transfer by mounting protrusions on surfaces 
The subject of enhanced heat transfer is of serious interest in heat exchanger applications. 
Concepts like space and energy optimization in power and process industries call for 
more and more compact designs of heat exchangers. Specially designed surfaces use 
protrusions mounted on them to enhance heat transfer. Some typical examples of the 
enhanced heat transfer surfaces which are very popular in different industrial applications 
(Webb, 1987) are shown in Fig. (2.2). There have been a number of recent survey articles 
and handbook sections (Bergles, 1978, 1983, 1985). Special surface geometries bring 
about the transport-enhancement by establishing a higher “hA” per unit base surface area. 
It may be mentioned that (1/hA) is the thermal resistance associated with heat transfer by 
convection at a surface. Some well-known methods of manipulation of surface-geometry 
are: 

1 . Inserting twisted tapes or turbulence promoters for internal flows in circular tubes. 

2. Roughening the surfaces for channel flows in closely spaced parallel plate channels, 
typical of plate-fin and plate tj^e heat exchangers. 

3. Improving flow structure in external flows normal to tubes and tube banks. 

2.1.3 Flows in Parallel Plate Channels Formed by Two Neighboring Fins of Plate- 
Fin or Fin-Tube Heat Exchangers 

It may be mentioned here that we are particularly interested in enhancement of heat 
transfer in the gas side of fin-tube heat exchangers (with flat fins) and enhancement of 
heat transfer in plate-fin heat exchangers. With this intent, we would like to study 
different investigations related to augmentation of heat transfer suitable for proceeding 
applications. 

Augmentation of heat transfer is of special interest in channel flows where the rate of the 
heat transfer between the fluid and the channel walls deteriorates as the boundary layer 
grows on the channel walls and the flow tends to become fully developed. Protrusions 
can be mounted on these channel walls in order to disrupt the growth of the boundary 
layer and thereby enhance the heat transfer between the flowing fluid and the channel 
walls. Two relevant applications using this kind of flow configuration are the heat trans 
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Pl= Axial distance through 



Figure 2.2 Various methods of making enhanced tubes: (a) helical rib on inner surface 
and integral fin on outer surface (b) internal fins on inner surface and porous coating on 
other surface (c) twisted tape insert on inner surface and integral fins on outer surface. 
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fer between the gas and the fin in the case of gas-liquid fin-tube cross-flow heat 
exchangers and the heat transfer between flowing fluid and plates in the case of plate-fin 
heat exchangers (Fig.2.3). The evolution towards a fully developed flow can be disturbed 



(d) (g) 


Figure 2.3 : (a)Plate -fin heat exchanger and its surface geometries: (b) plain rectangular 
fins (c) plain triangular fins (d)wavy fins (e) offset strip (f) perforated fins (g)louvered 
fins; after Webb (1987) 

by using a multi-louvered surface geometry for the plates. Investigations by Achaichia 
and Cowell (1988) provide a detailed performance data for louvered fin surfaces. 
However in using louvered fins, enhancement is obtained at the price of high-pressure 
drop. To circumvent this difficulty, protrusions in the form of slender delta wings or 
wingiets can be deployed (Fig.2.3). As shown, the base of the wing remains attached to 
the fin and the apex faces the incoming stream with an angle of attack with this 
configuration. The longitudinal vortices are generated along the side edge of the wing- 
shaped vortex generator due to the pressure difference between the front surface facing 
the flow and the back surface (Fig.2.3). These longitudinal vortices, generated by the 
vortex generators, can be made to disrupt the growth of boundary layer in a channel by 
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exchanging the fluid from the neax-wall-region with the channel-core-region and thus 
they can serve to enhance the heat transfer rate while producing less of a pressure drop. 
Use of longitudinal vortices for boundary layer control is well known (Pearcey, 1961) 
and the vortex generators are used in commercial airplanes for this purpose. 

MAC algorithm which is described below, calculates the velocitiy in U, V and W 
directions as well as it calculates the pressure P and temperature 9. These values are used 
by post-processing computer code (written separately), for the final calculation of Nusselt 
number. Apparent friction factor etc. 

2.2 Problem Formulation 

2.2.1 Introduction to the MAC (Marker and Cell) Method 

The important ideas on which the MAC algorithm is based are: 

1. Unsteady Navier-Stokes equations for incompressible flow in weak conservative 
form and the continuity equation are the governing equations. 

2. Description of the problem is elliptic in space and parabolic in time. Solution will be 
marched in the time direction. At each time step, a converged solution in space 
domain is obtained but this converged solution at any time step may not be the 
solution of the physical problem. 

3. If the problem is steady, in its physical sense, then after a finite number of steps in 
time direction, two consecutive time steps will show identical solutions. However in a 
machine-computation this is not possible hence very small value say, ‘STAT’ is 
predefined. Typically, STAT may be chosen between 10'^ and 10'^. If the maximum 
discrepancy of any of the velocity components for two consecutive time steps at any 
location over the entire space does not exceed ‘STAT’ tiien it can be said that the 
steady solution has been evolved. 

4. If the physical problem is basically unsteady m nature, the aforesaid maximum 
discrepancy of any dependent variable for two consecutive time steps will never be 
less than ‘STAT’. However, for such a situation, a velocity component can be stored 
over a long duration of time and plot of the velocity component against time depicts 
the character of flow. Such a flow may be labeled as simply ‘unsteady’. 
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5. With the help of momentum equations, we compute explicitly a provisional value of 
the velocity components for the next time step. 

Consider the weak conservative form of non-dimensional momentum equation in the X 
direction; 


du d{u^) ^ d(uv) ^ d{uw) _ dp ^ 
dt dx dy dz dx Re 

It is assumed that at t =nth level, we have a converged solution. Then for the next time 
step 


i,J,k 


At 


= [MSIDUl,j. 


or, 

[rESIDU]"jj^ consists of convective and diffusive terms, and the pressure gradient. 
Similarly, the provisional values for v~^ and w~^ can be explicitly computed. These 

explicitly advanced velocity components may not constitute a realistic flow field. A 
divergence free velocity field has to exist in order to describe a plausible incompressible 
flow situation. Now, with these provisional u~^ ,v~^ and values, continuity 

equation is evaluated in each cell. If V .V produces a nonzero value, there must be some 
amount of mass accumulation or annihilation in each cell, which is not physically 
possible, leading to an incorrect pressure field. Therefore, pressure at any cell is directly 
linked with the value of (V.V ) at the cell. Now, on one hand the pressure has to be 
corrected with the help of non-zero divergence value and on the other, the velocity 
components have to be adjusted. The correction procedure continues through an iterative 
cycle till the divergence free velocity field is ensured. 

Boundary conditions are to be applied after each explicit evaluation for the time step is 
accomplished. Since the governing equations are elliptic in space, boundary conditions 
on all confirming surfaces are required. Moreover, the boundary conditions are also to be 
applied after every pressure velocity iteration. The five principal kinds of boundary 
conditions to be considered are rigid no-slip walls, free-slip walls, in flow and outflow 
boundaries, and periodic (repeating) boundaries. (Biswas, G., (1995), Solution of Navi er- 
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Stokes Equations for Incompressible Flows Using MAC and SIMPLE Algorithms, in 
Computational Fluid Flow and Heat Transfer, ed. K. Muralidhar and T. Sundararajan, 
Narosa Publishing House, India.) 

2.2.2 Governing Equations for Laminar Flow 

The dimensionless equations for continuity, momentum and energy for laminar 
incompressible flow may be expressed in their conservative forms as following: 


Z) = 

dU dV 

dW 

= 0 



(2.1) 


+ 




dx'^ dY 

az 
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dU^ 

dUV 

dUfF 
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(2.3) 


+ + • 

+ 

“ 



dr 

dX 

dY 

az 

ar^ 

Re 

dW 

dUW 

dVW 

dW^ 

dP 

V^JV 

(2.4) 


-f 4" 


= 

+ 

dr 

dX 

dY 

az 

az 

Re 

ae 

due - 

dve 

dwe _ 



(2.5) 
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5t 

dX 

dY 

az 

Re.Pr 



In the above equations, velocities have been non-dimensionalized with the average 
incoming velocity Uav at the channel inlet, all lengths with the channel height H, the 
pressure with p and the non-dimensional temperature is defined as, 

Q = (X-T ^ )/ (Tw-T ^ ). The Reynolds numher and the Prandtl number are defined by Re 
and Pr respectively. 

2.2.3 Boundary Conditions for Laminar Flow 

The houndary conditions of interest in this investigation are (refer Fig. 2.4): 

• Top and bottom plates: 

U=V=W=0; e=\ (2.6) 

• Side wall (z = B/2) and midplane (z = 0): 
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fF= (— ) = (^) = o-(~) = 0 


(2.7) 


• At the channel inlet: 

U=U(YJ, V=fF=0;e = 0 ( 2 . 8 ) 

• At the exit; a smooth transition through the outflow boundary is ensured by setting: 

^ = -n 

aA" ax^ ’dX^~^ (2-9) 

No-shp boundary conditions for the velocities on the obstacles are used. The temperature 
of the obstacle is considered constant and equal to Tw 



k* kim 5x 


Out flow boundory 


Bollom boundory 


Right side woll 


Front inflow plane 


Figure2.4: Boundary conditions and fictitious boundary cells 
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2.2.4 Grid System Used 

Computational domain is divided into a set of rectangular cells (Fig.2.5(a)) and a 
staggered grid arrangement is used such that velocity components are defined at the 
center of the cell faces to which they are normal (Fig2.5 (b)). The pressure and 
temperature are defined at the center of the cell. In such an arrangement, pressure 
difference between two adjacent cells. The pressure field will accept a reasonable 
pressure distribution for a correct velocity field. Another important advantage of such a 
grid system is that transport rates across the faces of the control volumes can be 
computed without interpolation of velocity components. 

2.2.5 Boundary Conditions for Confining Walls 

Since, the equations are elliptic in space, we need boundary conditions on all confining 
walls, -even at the outlet. Consider, for example, the bottom boundary of the computing 
(physical) mesh. If this boundary is to be a rigid no-slip wall, the normal velocity on the 
wall must be zero. Here, we are considering a stationary wall. With reference to Fig.2.5 
(b), we have 




^.. 1 .* = 0 
^1,1, k ~ ~^i.2,k 


For i = 2 to ire 
And k = 2 to kre 


( 2 . 10 ) 


If the right side wall is a ffee-slip (vanishing shear) boundary, the normal velocity must 
be zero and the tangential velocities should have no normal gradient. 


=0 

U, J ^ = ~U, J 2 r For i = 2 to ire 
V =-V Andj=2tojre (2.11) 

' i,J,2 
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Figure 2.5 (a) Grid spacing in the computational domain and the location of the wmg type 
vortex grid showing the descretized variables (b) Three dimensional staggered grid 
showing the locations of the descretized variables. 
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The front plane is to be provided with inflow boundary condition. Any desired functional 
relationship may be recommended. Generally, normal velocity components are set to zero 
and a uniform or parabolic axial velocity may be deployed. Hence, with reference to Fig. 


(2.5(b)), we can write 


V =i-V 
Wuy.k = 

^.j,k = 1-0 

or,C/,^,,=1.5*[l-(^^)^] 


For j = 2 to ire 
y Andk = 2tokre (2.12) 


where, jm is the horizontal midplane. 

Continuative or outflow boundaries always pose a problem for low-speed calculations, 
because whatever prescription is chosen it can affect the entire flow field upstream. What 
is needed is a prescription that permits fluid to flow out of the mesh with a minimum of 
upstream influence. In MAC family of algorithms, the second derivatives of the 
dependent variables in flow direction are set to zero m order to ensure smooth transition 
through the outflow boundary. These are implemented in the following way; 


Ui+ij^k Ui-ijk 

Vi+ij_k = 2Viik—Vi.ij^k ^ 
Wi+ij_k = 2Wij^k— Wi-ij^k 


For j = 2 to jre 
And k = 2 to kre 
At i = lim - 2 


(2.13) 


2.2.6 Boundary Conditions for the Obstacle 


A special effort is required to satisfy the kinematic boundary conditions on the vortex 
generator. From Fig. (2.6), the angle of attack of the wing-type vortex generator can be 
expressed as j3 = tan'^ {5Y ! 5X) . Fig. (2.6) further reveals that the U and V components 
of velocity fall directly on the surface of the obstacle in the region of the wing and are 
therefore equal to zero. Implementation of no-slip condition for U component velocity at 
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the side edge of the obstacle needs some manipulations. From Fig.(2.6), it is found that 
Uij,k+i /2 = 0. In the other words, we can write 

Uij,k“ "Uij.k+l 

(for i=ia to ib, j =2 to j a and k = 2 to kc) (2.14) 

It may be mentioned that the difference of count between kc and 2 is equal to that 
between ib and ia. The W component of velocity also needs manipulation for its zero 
value on the surface of the obstacle. From Fig.(2.7) 



Figure 2 6: Relative location of velocity components and wing on X-Y plane at Z=0 



(for i = ia to ib, j =2 to ja and k =2 to kc) 

The temperature boundary conditions on the wing-type obstacle are diagrammatically 
shown in Fig.(2.8). The obstacle boundary conditions for temperature are: 
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Figure 2.8: Thermal boundary conditions on the wing-type vortex 










(for i =ia to ib, j = 2 to ja and k =2 to kc) 

In practice, the vortex generators can easily be manufactured by punching or embossing 
the wall. In some computational results, the effects due to the hole beneath the vortex 
generator have been taken into account. It is considered that the heat exchanger cores 



Figure 2.9 : Periodic boimdary conditions for pressure and velocities on the top and the 
bottom wall. 
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have a number of plates and on each plate the vortex generators are punched out in such a 
way that the holes are perfectly aligned. In our computational domain, a small part of top 
and bottom boundaries must be set to identify the punched holes. For this purpose 
spacewise-periodic boundary conditions on the location of the punched holes have to be 
applied in the y-direction (Fig.2.9). Here, the period length is channel height ‘if. 
However boundary conditions for the fictitious cells on the bottom are; 


and on the top 





~~ ^i,jre,k 


= 0.5(V,jre.k+V,M) 


= Kjre.li 


= P 

ijre.k 




(for all i, k in the punched holes) 


(2.17) 


^ i, Jim, k 


y 


V 

ijre,k 

= 0.5(F,, 

w 

i,jm,k 

= 

P 

i,jm,k 

~ ^i,2,k 


(for all i,k in the punched holes) 




y 

J 


(2.18) 


Such boundary conditions have also been discussed by Hirt et al. (1975) and Biswas and 
Chattopadhyay (1992). However, for plate-fin heat exchangers, a different manufacturing 
technique for mounting the vortex generators is adapted. For this case, no hale exists 
underneath the vortex generators in order to avoid mixing of the hot and cold streams. 
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2.3 MAC (Marker and Cell) Algorithm 

Because of staggered grid arrangement (Fig.2.5), the velocities are not defined at the 
nodal points. Whenever required, the velocities at the nodal points are to be found by 
interpolation. For example, with uniform grids, we can write Ui- = 0.5(LJi.ij,k+Uij,k). 
Where a product or square of such quantity appears, it is to be averaged first and then the 
product is to be formed. 

Convective terms of the momentum equations are discretized using a weighted average of 
second upwind and space centered scheme (Hirt et al., 1975). Diffusive terms are 
discretized by a central difference scheme. Let us consider the discretized terms of x- 
momentum equation (refer to Fig. *). 


'dU^~ 

1 

‘(V,,,. + V..,,,. XV.,., + ) + « , 1 (V.,.. + ) 1 (V J,* - V. w ) ■ 

- . 

i,j,k 



= DUUDX (2.19) 


'dUV' 


1 

+v.,.u)+“, -V j.u) 

. sx . 

ij.k 

~ 45Y 



=DUVDY (2.20) 


'dUW' 

1 


. dZ _ 




^DUWDZ (2.21) 


dx ' 

d^U 

dX^ 


i-^\,j,k ^ i,j,k 


= DPDX 


5X 

(5A^) 


= D2UDX2 


( 2 . 22 ) 

(2.23) 
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with, a p > 1, the scheme tends to second upwind and > 0, the scheme tends to 

Space centered. 

The factor is chosen in such a way that the differencing scheme retains ‘something’ of 

second order accuracy and the required upwinding is done for the sake of stability. A 
typical value of is between 0.2. and 0.3. 

Solutions for the velocities are obtained in two folds. First, the velocity components are 
advanced explicitly using the previous state of the flow, having calculated accelerations 
caused by convection, diffusion and pressure gradients through a time step St . The 
explicit time increment may not necessarily lead to a velocity field with zero mass 
divergence in each cell. In the subsequent fold, adjustment of pressure and velocity is 
done by an iterative process in order to ensure mass conservation in each cell. If the 

explicitly advanced provisional velocity is termed as U”^ , then from equation (*) we 
can write 

+ St[XESIDUI,. (2.24) 

where. 


[RESIDUl:,,, = 


i-DUUDX - DUVDY - DUWDZ) - DPDX 
+ (1 / P£){D1UDX2 + D2UDY2 + D2UDZ2) 


In a similar manner we evaluate 


(2.25) 




(2.26) 


K'* ’Kj: +H«ESIDWI,j, (2.27) 

As discussed earlier, the explicitly advanced tilde velocities may not necessarily lead to a 
flow field with zero mass divergence in each cell. This implies that at this stage the 
pressure distribution is not correct. Pressure in each cell will be corrected in such a way 
that there is no net mass flow in or out of the cell. In the original MAC metohd, the 
corrected pressures are obtained from the solution of a Poisson equation for pressure. A 
related technique developed by Chorin (1967) involved a simultaneous iteration on 
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pressures and velocity components. Vicelli (1971) has shown that the two methods as 
applied to the MAC method are equivalent. We shall make use of the iterative correction 
procedure of Chorin (1967) in order to obtain a divergence-free velocity field. The 
mathematical methodology of this iterative pressure- velocity correction procedure will be 
discussed herein. 

The relationship between the explicitly advanced'velocity component and the velocity at 
the previous time step may be written as 

~ Slip” + P” 1 

— rrn . t i,j,k rl'< '^o^ 


rjn+\ _rjn 


■ + 5x[CONDIFUXj 


[CONDIFUXjj^ = [RFSID + [dpdx 


(2.28) 


(2.29) 


where, on the other hand, the corrected velocity-component (unknown) will be related to 
the correct pressure (also known) in the following way: 

•"■'■‘‘ +5t[CONDIFUI^_, (2.30) 


■ + 5r[CONDIFU]ljj^ 


(2.30) 


from equations (*) & (*) we can write 

rrn+l frn+l _ ~ '^+ 1 . 7 .* 1 


i,j,k ^ i,j,k 


(2.31) 


where the pressure is defined as 


p' -- p"+l _ p” 


(2.32) 


Neither the pressure corrections nor are known explicitly at this stage so one can 
not be calculated with the help of the other. Calculations are done in an iterative cycle 


and we write 




\pi,J,k ^i*l,j,k ] 


In a similar way, we can formulate the following array: 


Corrected Estimated 


Correction 


^AKj,k 


(2.33) 
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The correction is done through continuity equation. Plugging-in the above relationship 
into the continuity equation (*) we shall obtain 
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+ 




5X^ 


- + 


§7" 


5Z' 


(2.40) 


In deriving the above expression, it is assumed that the pressure corrections in the 
neighboring cells are zero. Back to the calculations, we can write 




05^/1 1 1 

25x (— — :r + 1 — :r + 
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5X^ 5Y^ 5Z\ 


os / 1 1 1 
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(2.42) 


In order to accelerate the calculation, the pressure correction equation is modified as 


P = — 

^ij,k 




26t( 
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1 


- + - 


1 


5X^ 5Y^ 5Z^ 


(2.43) 


where ©ois the overrelaxation factor. A value of © 0 = 1.7 is commonly used. After 


calculating P,'^ the pressure in the cell (i,j,k) is adjusted as 

P.7!k^P:j.k+Kj.k (2-44) 

Now the pressure and velocity components for each cell are corrected through an iterative 
procedure in such a way that for the final pressure field, the velocity divergence in each 
cell vanishes. The process is continued till divergence-free velocity field is reached with a 
prescribed upper bound; here a value of 0.0001 is recommended. This solution scheme is 
continued until a steady flow is obtained. 

Finally we discuss another important observation. If the velocity boundary conditions are 
correct pressure will be determined in all the cells including the cells at the boundary. 
Thus, this method avoids the application of pressure boundary conditions. This typical 
feature of modified MAC method has been discussed in more details by Peyret and 
Taylor(1983). 
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Discretization Technique for k and s Equations: 

Convective terms of k and s equations are also discretized using a weighted average of 
second upwind and space centered scheme (Hirt et al., 1975). Diffusive terms are 
discretized by a central difference scheme. 

Let us consider the discretized terms of the k equation 


~dkL/' 
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. 3X _ 

ij.k 

" 25Y 

i ^z-l,y,yfc 1 (^z-l,y,it ~’^z,y,)t)_ 
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^DIFFkX (2.46) 

The difference quotients involved in the Generation term are discretized as following. 
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( 2 . 51 ) 
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dJV 1 


dX 4SX 
= DWDX 


(2.52) 




= DUDZ 


(2.53) 


5F 1 


dZ 45Z 
= DVDZ 


k...+Vu+’', ,7,^+1 + ^..V 


(2.54) 


dJV 1 




dr 45X ' 

= DWDY (2.55) 

In terms of the difference quotients, the s equation consists of the similar quotients as 
the k equations. 


Solutions of the Energy Equation: 

After evaluating the correct velocities, the energy equation is solved with an successive 
Over-Relaxation technique to determine the temperature field. We shall discuss the 
solution procedure of the energy equation in this section. 

The steady state energy equation, neglecting the dissipation term, may be written in the 
following conservative form as 


due dve dwe _ v^e 

dX dY dZ ~ Pe 
Equation (2.56) may be wntten as 

V^e = Pe[CONVTXj^^ 


(2.56) 

(2.57) 
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where, \C0NVT\ ^ ^ is the discretized convective terms on the left hand side of equation 

(2.56) and m stands for iterative counter. To start with, we can assume any guess value of 
6 throughout the flow field. U, V and W are known from the solution of momentum 
equation hence equation (2.57) is now linear equation. However, from the guess values of 
6 and known correct values of U, V and W the left hand side of equation (2.56) is 
evaluated. A weighted average scheme may be adapted for discretization of the 
convective terms. After discretizing and evaluating right hand side of equation (2.57) we 
obtain a Poisson equation for temperature with a source term on the right hand side. Now, 
we can follow a SOR technique for solving equation (2.57). Consider a discretized 
equation as 


, ^ij+l,k i,j,k ^ t,j-l,k , ^ tj.k ^ t,j,k-l _ p»m 

■ "i 1 — u 


i5xy 
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where. 


(5J7)" (57)" (5Z)" 

0, ^ 4 in equation (*) may be assumed to be most recent value and it may be written as 
0," . In order to accelerate the speed of computation we introduce a over-relaxation 
factor CO and 


0 m+l 
hJ,k 






J,k 



(2.59) 
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where the previous value, 0"^ ^. is the most recent value and 0”*| is the calculated 

better guess. The procedure will continue till the required convergence is achieved. 


Numencal Stability Considerations: 

For accuracy, the mesh size must be chosen small enough to resolve the expected spatial 
variations in all dependent variables. 

Once a mesh has been chosen, the choice of the time increment is governed by two 
restrictions, namely, the Courant-Friedrichs-Lewy (CFL) conditionand the restriction on 
the basis of grid-Fourier numbers. 

According to the CFL condition, material cannot move through more than one cell in one 
time step. Therefore, the time increment must satisfy the inequality 


5x <min' 


5X 5Y 

\U\’\VY\W\] 


( 2 . 60 ) 


where the minimum is with respect to every cell in the mesh. Typically, 8x is chosen 
equal to one-fourth to one-third of the minimu m cell transit time. When the viscous 
diffiision terms are more important, the condition necessary to ensure stability is dictated 
by the restriction on the grid Fourier numbers, which results in 


^ 1 {5Xy{8Y)\5Zy 
2{8X^ +8Y^ +8Z^) 


(2.61) 


The final 8x for each time increment is the minimum of the 8x ’s obtained from 
equations (2.60) and (2.61). 

(This part of the thesis has been taken from Deb, Prashanta, (April, 1994), “Three 
Dimensional Computation of Laminar and Turbulent Flows with Heat Transfer in a 
channel with Built-in Longitudinal Vortex Generators” , PhD. Thesis, ME Department, 
IIT Kanpur.) 


2.4 Literature Survey 

What originally motivated the introduction of EVOP, was the idea that the widespread 
and daily use of simple statistical design and analysis during routine production by 
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process operatives themselves could reap enormous additional rewards. [Box and Draper 
(1969)] 

An excellent review of industrial applications of EVOP has been presented by Hunter and 
Kittrell (1966). This review, which lists 68 references, discusses applications of EVOP in 
a wide variety of environments. Among these, applications in the chemical industry are 
most numerous, with special references to uses by the Dow Chemical Company 
[Anonymous (1961u, 1961^>)], Imperial Chemical Industries Limited [Coutie (1959a, 
1959b)], Tennessee Eastman Company [DeBusk (1962); Pursglove (1961); Wilson 
(I960)], the Chemstrand Corporation [Annual Report (1959); Riordan (1958, 1956b)]’ 
Phillips Petroleum Company [Pmsglove (1961); Weaver (1963)] 

, Standard Oil of Ohio [Klingel and McIntyre (1962); Pappas(1962)], and Monsanto 
Company [Anonymous (1960); Hehner (1963c)]. Several applications in the food 
industry are also mentioned, in particular use by Swift and Company [Samuel (1962)], 
Canadian Packers Limited [Samuel(1962)], and the A. E. Stanley Manufacturing 
Company [Koleff (1963)]. Also discussed are uses of EVOP by the canning industry 
[Filice (1963)], and by the Maumee Chemical Company in such diverse projects as 
saccharin production, biocide for sea lampreys isatoic anhydride, anthranilic acid and 
benzotriazole processes [Anonymous (1963)]. Other applications are mentioned in the 
automotive industry where EVOP was applied to resistance welding of automotive sheet 
metal [Thomas and Webster (I960)] and other manufacturing operations [Thomas 
(1965)]. Among these applications a savings of several hundreds thousand dollars per 
year by the Chemstrand Corporation has been mentioned by F. S. Riordan (1958,1959b), 
and a savings of $250,000 per year for two products has been reported by H. O. Hehner 
(1963c) in the Monsanto Company. In the description of uses by R. E. DeBusk (1962) 
that EVOP was used in approximately 15 processes with a savings of from $15,000 to 
$50,000 per year per process. 

Box and Draper (1969) have regarded EVOP as a simple but powerful statistical tool with 
wide application in industry. Their experience has long shown that statistical methods, 
sometimes quite sophisticated in character, can be of great value in improving the 
efficiency of laboratory and pilot-plant investigations made by specially trained chemists 
and engineers. 
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2.5 Scope of the Present Work 

Since the details of the interactions between the ‘angle of attack and length of winglef 
with ‘enhancement of heat transfer and loss of pressure energy’ are not still clearly 
understood, there is a general agreement that more research on this subject is required. 

In this present work, a detailed three-dimensional numerical model of the MAC 
Algorithm is used which is able to provide a clear understanding of the physics of the 
flow and heat transfer. With the numerical model full Navier-Stokes equations together 
with governing equations for energy have been solved in a rectangular channel with built- 
in winglet vortex generators. The present study is aimed at a numerical investigation of 
heat transfer enhancement in a channel with built-in winglet type vortex generators using 
EVOP technique combined with Multi-criteria Decision Making and Pareto Optimality 
concepts. 

The objectives of study can be explicitly stated as below: 

• To optimize the delta winglet geometry in the channel so as to enhance the heat 
transfer by 30 %, or more, than that of the plain channel. 

• To optimize the pressure drop m the channel flows and keep it as close as possible to 
that of the plain channel flow. 

For the present study, Reynolds number is 790 and Prandtl number is 0.7, for all the heat 
exchanger design evaluations. 
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3. Heat Exchanger Optimization Formulation 


3.1 Problem Optimization by EVOP 

For the optimization of the problem above, there are two input variables and two output 
variables. The input variables are the 'angle of attack" and 'the length of the winglef 
(refer Figure 3.1 and Fig. 3.2). 



Figure 3.1 : Geometry of Channel 



where, 

B= width of channel, P = angle of attack, H= height of channel, 1 = length of winglet, 
L= length of channel, b = height of winglet 


39 


The two performance criteria that we wanted to optimize are ‘Nusselt number’ and 
‘Apparent friction factor’. For the EVOP technique, we calculated the average Nusselt 
number spanwise and further averaged it using ‘Simpson’s one third rule’ of integration 
(Fig. 3.3). Thus, unique Nusselt number for each angle of attack its respective length 

of wmglet was calculated. This averaged Nusselt number was used as the measure of the 
heat transfer for the channel. Pressure drop was also calculated using average method. 
The pressure at the entry of the flow was calculated in the ‘j ’ and ‘k’ directions and is 
averaged. The pressure at exit was also calculated and averaged. The difference between 
these two spanwise average pressures gives us the ‘pressure drop’. 

This ‘pressure drop’ was chosen as the second criterion for the optimization. So the task 
before us now was to optimize the input design variables, "angle of attack" and "length of 
winglef in such a way that it would give us maximum Nusselt number while producing 
mimmum pressure drop. 

The EVOP discussed in the Chapter 1 can in general deal with ‘n’ number of input 
variables, but it takes into account the single output parameter, for optimization. The 
special feature of the present study is that, by varying the two input variables, vise "angle 
of attack" and "length of winglef, we maximize the response Nusselt Number while at the 
same time we are trying to keep the pressure drop as close as possible to that of a plain 
channel. To achieve this we employ ‘Pareto Optimality’ (Section 1.5) concepts. 
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Figures. 3 (a): Symbolic Spanwise Average Nusselt Number Distribution 



Figure 3.3 (b): Heat Exchanger with Upper and Bottom Plates showing cells in flow 

direction. 


* iim is the total number of cells in the direction of flow. The first and last cell in 
the geometry are fictitious. 
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3.2 Methodology for Heat Exchanger Geometry Optimization 


Let us consider Figure 3.4 as shown below. Here the Nusselt number is plotted on the X- 
axis, and the Pressure Drop is plotted on the Y-axis. The task ahead of us is to maximize 
Nusselt number, of the channel with winglet, while keeping the Pressure Drop as close as 
possible to that of the plain channel. 

Let the points shown on the graph represent the outputs of the different runs for various 
‘angle of attacks and length of winglet’ (design variables) settings. Then from the figure 
we see that the points numbered 1,2, and 3 are representing the first ‘Pareto front’. So, the 
neighborhood of these angles of attack and the length of winglet values are to be 
searched. 



Nusselt Number 

Figure3.4; Symbolic representation of Bi-criteria problem solutions in the Response 

Space 

Suppose we take point ‘1’ to start the EVOP study . Refer to Figure (3.5) below, which 
shows the four variants of point ‘1’ namely la, lb, Ic and Id. In this figure, the angle of 
attack is shown on ‘X’ axis and the length of winglet is shown on ‘Y’ axis. 

The same EVOP procedure is to be carried out for points 2 and 3. All new points are to 
be tested for the Pareto optimality. So, with the new inputs the process is run. These new 
inputs will give us new output values. The new output values are now again plotted on 
the graph to identify the members with Pareto property. (Figure 3.5) 

The same procedure will be carried out for points 2 and 3 also. The new points obtained 
are tested for the optimality. So, with new inputs the process is ealuated. These new 
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inputs will give us new output values. The new output values are now again plotted on 
the bi-criteria response space. 
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Figures. 5: Two dimensional variants in the Design Space 


Refer to Figure (3.6) shown below. 

By looking at the figure we see that a new (second generation) Pareto front consisting of 
Id, 2a, 2b and 3a has emerged. 



Nusselt Number. 

Figures. 6: New Pareto front in the Response Space after new input runs 
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The new Pareto front consists of points Id, 2a, 2b and 3a. These points are plotted again 
on the ‘angle of attack versus length of winglet’ graph, and new EVOP variants are 
selected again and the process is re-run. 

This process of plotting output and then selecting new input parameters accordingly is 
conducted iteratively, till the Pareto front is saturated and does not improve much. After 
scituration is reached, the senes of input variable pairs representing the solution on the 
Pareto front will give us the set of optimal solutions, or input parameters. 

We must point out here that the process of optimization involves the judgement of the 
designer. Thus the decision to fathom out some solution area or exclude some input from 
further searching is really subjective. So this process of EVOP-based design optimization 
would depend on experience, knowledge and the ingenuity of the designer himself 

3.3 Determination of Performance Parameters 

Calculation of Optimization Parameters: 

The Spanwise Average Nusselt number is calculated. Fig. 3.3(a) shows the symbohc 
curve for the Nusselt number distnbution spanwise, along the direction of the flow i.e. I 
direction. To calculate the spanwise average of Nusselt number the following procedure 
is adopted: 

• The Nusselt number at plane number 1 is calculated on bottom plate and is averaged 
on the ‘K’ direction. 

• The Nusselt number at plane number 1 is calculated on upper plate and, is also 
averaged on the ‘K’ direction. 

• Now, the average of these, bottom plate Nusselt number and the upper plate Nusselt 
■ number is calculated. This is spanwise average Nusselt Number. 

Fig.3.3(a) shows the spanwise Nusselt number, at position ‘1’, shown in Fig.3.3(b). 

The area under curve shown in Fig. 3.3(a) is calculated using, Simpson’s one-third rule of 
integration. The final Nusselt number, which is used as one of the optimizing objective, is 
calculated by averaging this area on the length of the channel in the direction of flow. 

The pressure drop is calculated as follows: 
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• Spanwise average apparent friction is calculated at the entry of the heat exchanger, 
along the plane ‘1’, where the node number is 2 in the flow direction. 

• Spanwise average friction factor is calculated at the (iim-1)^ plane. 

• The subtraction of apparent friction factor at plane 1, and the apparent friction factor 
at plane (iim-1) will give us the value of the second objective of optimization. 


' Where, iim iVthe total number of cells m the direction of the flow m the channel. The first and last cells 
are fictitious 
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4. Results and Analysis 


4.1 Initial Solutions 

4.1.1 Channel Geometry 

EVOP technique specifically operates on improving process after its start-up. So to get started 
with initial solutions we must first of all fix up the geometry of the channel. Subsequently 
optimization will be carried out in a direction to enhance the Average Nusselt Number while 
maintaining the pressure drop to its minimum possible. 

In selecting the Channel Geometry the dimensions are so chosen that it will accommodate a 
change in length of winglet up to 1 50%and a change in angle of attack from 0° to 90°. 



(c) 


Figure (4.1): Channel Geometry: (a) Geometry of heat exchanger plates (b) Geometry of 
winglet (c) Top view of the winglet on heat-exchanger plate 
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The following dimensions were used in the numerical investigations of this study; 

B = 4.0, L = 1 1 .5, H = 1 .0, Li = 2.5 (vanable), L 2 = 4.0, Bi = 3.0 
The purpose behind choosing this kind of geometry was that by fixing the trailing edge position 
of the winglet, any changes in the angle of attack and the length of winglet could be 
accommodated. This will lead to draw a fair conclusion about the effects of inputs on the 
outputs, and the decision making task could be eased. 

As shown in Figure 4.1, all the dimensions related to Channel and the position of the trailing 
edge of the winglet are fixed. 

4.1.2 Sample Calculations for determining Input Parameters for MAC ( Marker and Cell) 
Algorithm 

The inputs to the MAC algorithm for numerical run are calculated as shown below. For 
demonstrating the sample calculation the angle of attack is chosen equal to 30° and the length of 
winglet is chosen equal to 2.5. 



Figure 4.2; (a) Channel dimensions 
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Figure (4.2): (b) Sample Winglet Dimensions 

After fixation of the input parameters, the calculation of input parameters for MAC algorithm are 
carried out as shown below: 

fy = 1.0, Ay = 0.05, p = 30°, H = 1.0, B = 4.0, L = 11.5, L 2 = 4.0, Bi = 3.0. 


Az = 4*0.05 
fj = Li * sin p 

.•.f, = 2.5*sin30 = 1.25. (4.1) 

from' equation (4.1) 

Az = 1 .25 * 0.05 = 0.0625 (4.2) 

(jim-2) * Ay = H 

.-.jim = 22 (4-3) 

P = tan"' (Az / Ax) 

.-. Ax = (0.0625 / tan 30) = 0.108 (4.4) 


from equation (4.2) 

B = (kim-2) * Az = 4.0 

.■.kim = 66 (4.5) 

from equation (4.4) 
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L = (iim-2) * Ax 
.•.iim =108 

kc = (Bi / Az) + 1 
.•.kc = 29 


ib = (Li / Ax) + 1 
.■.ib=37 


(4.6) 


(4.7) 

+ 1 

(4.8) 


(4.9) 


kb = {(Bi - Li sin p) / (Az)} 
.•.kb=9 


H= {(ib-ia)+l} * Ay =1.0 


.•.(ib-ia) = 19 


(4.10) 


from equation (4.9) and (4.10) we calculate, 
.•.ia=19 (4.11) 


where, 

iim = total number of cells in ‘i’ direction, (flow) 
jim = total number of cells in ‘j’ direction, (height) 
kim = total number of cells in ‘k’ direction, (span) 
kc = coordinate of trailing edge of winglet in ‘k’ direction, 
kb = coordinate of leading edge of winglet in ‘k’ direction, 
ib = coordinate of trailing edge of winglet in ‘i’ direction, 
ia = coordinate of leading edge of winglet in ‘i’ direction. 

Finally, the values of iim, jim, kim, fx, fy, fz, kb, ia and ib are the inputs to MAC algorithm. The 
procedure elaborated above, is carried out for every new input angle of attack and length of 
winglet. 
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L = (iim-2) * Ax 
.'.iim =108 


(4.6) 


kc = (Bi / Az) + 1 

.•.kc = 29 (4 7 ) 

kb = {(Bi - Li sin P) / (Az)} + 1 

.•.kb=9 (4.8) 

ib = (L| / Ax) + 1 

•■•ib=37 (4.9) 

H = {(ib - ia)+l} * Ay =1.0 

.•.(ib-ia) = 19 (4.10) 

from equation (4.9) and (4.10) we calculate, 
.•.ia=19 (4.11) 


where, 

iim = total number of cells in ‘i’ direction, (flow) 
jim = total number of cells in ‘j’ direction, (height) 
kim = total number of cells in ‘k’ direction, (span) 
kc = coordinate of trailing edge of winglet in ‘k’ direction, 
kb = coordinate of leading edge of winglet in ‘k’ direction, 
ib = coordinate of trailing edge of winglet in ‘i’ direction, 
ia = coordinate of leading edge of winglet in ‘i’ direction. 

Finally, the values of iim, jim, kim, fic, fy, fz, kb, ia and ib are the inputs to MAC algorithm. The 
procedure elaborated above, is carried out for every new input angle of attack and length of 
winglet. 


49 



4.1.3 Calculation of Output (Response) Parameters 


EVOP, would require responses, Nusselt number and pressure drop to be single-valued outputs 
for a particular angle of attack and length of winglet given as inputs. This was determined as 
follows. 

For the calculation of pressure drop the spanwise averaged apparent friction factor was used. 
The difference between the spanwise averaged apparent friction factor at flow entry and that at 
flow exit was calculated which provided us the pressure difference between entry and exit. This 
is called the “pressure drop” throughout the current optimization process. For comparison, the 
pressure dro without any winglet is 0.01 (in the apparent friction factor scale). 

For the calculation of average Nusselt number, first from MAC algorithm spanwise average 
Nusselt number at every cell along the flow direction was calculated. This gives us the spanwise 
distribution of average Nusselt number. In order to reach to a single unique number, Simpson’s 
one-third rule of integration was used. (The detailed discussion of calculation of input and 
output parameters is done in Chapter 3.) In order to achieve this, separate code was written and 
was included in the post-processing code for MAC algorithm. 

4.1.4 First phase of EVOP 

To get started with the length of winglet, Li was chosen to be equal to 2.5. We treated the MAC 
algorithm as a ‘Batch Process’, in which two process variables or factors could be manipulated. 
These two factors are the angle of attack and the length of winglet. In order to generate a set of 


0 10 20 30 40 

Angle of Attack 

Figure 4.3: EVOP Table for Initial Runs in 
the input space 

solutions, the length of winglet was fixed at 2.5, and by changing the angle of attack wise 15", 
22.5°, 30°, 37.5° the numerical runs were carried out.(We note that each of these MAC runs took 


"Hb 

a 


4h 

o 


0) 


3 

2 .5 
2 

1 5 
1 


50 



over 3 hours on the SGI Origin 2000 system.) The EVOP table for this initial set of inputs is 
shown in Figure 4.3. 

After the initial set of numerical mns were conducted, we got the two responses as shown in 
Figure 4.4, in which Nusselt Number is plotted on X-axis and the Pressure Drop is plotted on the 
Y- axis. This gives us the initial set of solutions. These solutions were treated as the starting 
points for the EVOP optimization. 



Figure 4.4: Initial Performance 


Based on these outputs we have to take the decision regarding the next EVOP runs, in order to 
go in a direction of optimizing them. 

4.2 Decision Making with EVOP 

From the ‘Information Board’ (i.e. Fig. (4.4)) we see an useful pattern emerging from the 
increment in angle of attack while the length of winglet was held constant. The increment in 
angle of attack improved the Average Nusselt Number, but at die same time, it also increased the 
pressure drop. So, to get a picture of what a change in length of winglet brought about, we ran 
the MAC algori thm again for same angles of attack in the first phase of EVOP. But this time, the 






length of winglet was increased up to its maximum feasible length in this particular chosen 
geometry of channel.' 

By changing the length from Li = 2.5, to Li = 3.75 in the second phase of EVOP, the numerical 
runs were carried out at the angles of attack 15°, 22.5°, 30° and 37.5° respectively. Fig. 4.5 (a) 



0 10 20 30 40 

Angle of Attack 


Figure 4.5 (a) EVOP table for second phase 

shows the EVOP table for second phase. And Fig. 4.5(b) shows the result table for after the 
numerical run with second phase of EVOP. 

From the ‘Result Table’ in Fig. 4.5(b) we infer that, increase in the length of winglet also has the 
same impact on the pressure drop. And also we can see that though liie pressure drop has 
increased over the first phase runs, the Nusselt Number also has increased by a considerable 
mar gin So, from this behavior we can conclude that by increasing angle of attack or the length 
of winglet we are going to get higher Average Nusselt Number accompamed by increased 
Pressure Drop. 

' The channel geometry is designed in such a maimer that, it can accomodate 150% change m length of winglet, that 
was used in the initial runs. 
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Next we explore the effect of reduction in the length of winglet for the set of angle of attack 



Figure 4.5 : (b) Results after second phase of EVOP 

in initial set of numerical runs. We would expect (before the actual numerical runs are carried 
out), that the pressure drop will be lower and the average Nusselt number also will be lower than 

for 

the set of inputs in the initial runs, if we decrease the length below 2.5. As we have fathomed out 



Figure 4.6; (a) Table for 3 phase of EVOP 

the upper limit of the length of winglet. we shall now also fathom out the lower limit of the 
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same. So arbitrarily we chose the length of wiglet to be equal to 1.25 (i.e. 50% of that of the 
initial inputs). As from the previous two phases of EVOP we can see that 15° angle of attack is 
not giving us any promising results, we eliminated it in this run, and the input angles of attack 
that we chose for third phase were 30° and 37.5°. The new EVOP table is as shown in Fig. 4.6(a). 
And the Result Table is shown in Fig. 4.6(b). From the Result Table we see that, the pressure 
drop for these new runs (marked by triangles) has decreased drastically and also the Nusselt 
number has decreased. Both the points are dominated by point ‘1’. 

0.02 
0.0B 
0.0B 
0.017 

I 0.016 

I 

I 0.014 
0.013 
0.012 
0.011 
0.01 


Figure 4.6: (b) Results after 3'"'^ phase of EVOP 
The Results after third phase are summarized in the Table 4.1. 
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4.3 Analysis of Solutions 


As the three phases are over we now have established the relationship among input parameters, 
angle of attack and the length of winglet with the output parameters wise Average Nusselt 
number and the Pressure drop. At this stage we have to decide on the further input sets in order 
to go in a direction towards the objectives of EVOP optimization. To take the decision about it, 
Pareto Diagram was plotted as shown in Fig. 4.7, which shows the Pareto Diagram with 
separately indicating the set of Pareto solutions so far. These solutions are picked out for the 
further decision of new input sets. 


Angle of 
Attack 

Length of Winglet 

Nusselt Number 

Pressure Drop 

15 

2.5 

13.37 

0.0103 

22.5 

2.5 

14.66 

0.015 

30 

2.5 

15.35 

0.016 

37.5 

2.5 

15.68 

0.017 

15 

3.75 

14.31 1 

0.013 

22.5 

3.75 

15.45 

0.019 

30 

3.75 

16.76 

0.035 

37.5 

3.75 

17.34 

0.038 

30 

1.25 

12.34 

. 0.011 

37.5 

1.25 

12.49 

0.0125 


Table 4.1 : Summary of results after 3’^'* phase of EVOP 
Table 4.2 enlists the Pareto solutions shown in Fig. 4.7. 

After studying the graph shown in Fig. (4.7), we have decided not to fathom out the angles of 
attack wise, 15° and 22.5°. The reason for not fathoming 15° anymore is, for the largest length of 
winglet also, the Nusselt number is not improving much. And for the other lengths of wmglet, 
the Nusselt number for 15° angle of attack, will lie somewhere below the 15° and 3.75 winglet 
length point on ftie result table. The reason for discarding 22.5° angle of attack is that, the output 
of 22.5° angle of attack with winglet length 3.75, is dominated by other points. And there is no 
need to fathom out this angle of attack further, as the Nusselt number is not going to improve 

further as well. 
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Pressure Drop 
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Figure 4.7; Pareto Solutions 



Table 4.2: List of Pareto Solutions after 3^“* phase 


We plotted EVOP table for these solutions. (Refer Fig. 4.8(a) to Fig. 4.8(e)) These EVOP tables 
have now evolved the further input sets which we have to fathom out m the 4^'' phase of EVOP. 
Table 4.3 enlists the new set of design inputs for the 4* phase. The decision of substitution of 
older input variables with the new input variables has been taken intuitively. 
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This new set of design inputs was fathomed out for its Pareto optimality by running the MAC 
algorithm with these inputs. 

The procedure, that is explained above, for the determination of the new set of input variables 
was carried out iteratively. The final set of Pareto Solutions is to be fathomed further. 




Variation introduced for 30° angle of Variation introduced for 37.5° angle of 
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Variation introduced for 37.5° Angle of 
Attack 


(d) 

Variation introduced for 30° angle of 
Attack 


Figure 4.8; EVOP Table for determination of Input sets of 4’^*' phase 
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Angle of Attack 

Length of Winglet 

40 

2.5 

37.5 

3.125 

37!5 

3.5 

40 

3.75 

30 

2.25 

30 

3.125 

30 

3.5 


Table 4.3: Set of new inputs for 4^^ phase of EVOP 

4.1 Results of EVOP Optimization 


EVOP optimization was run iteratively and the final set of optimal solutions was reaached. In 
deciding the optimality, the lower bound on the Nusselt number was considered to be 15. 
According to this condition the optimal set of solutions was marked. Figure 4.9 shows the output 
of all the runs that were run for the optimization process, and Figure 4.10 shows the EVOP 
Diagram for the entire optimization process. The Pareto optimal solutions are also marked in 
fig.(4.11)andfig.(4.10). 



Nusselt Number 


♦ length = 2.5 

♦ length = 3.75 
A length = 1 .25 
-length = 3.125 
X length = 3.5 
A length = 3.25 

♦ length = 2 25 


Figure 4.9: Output of all Evaluations 
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Figure 4.10: Pareto Optimaland Non-Pareto Design outputs 




Figure 4.1 1: EVOP Table for Pareto Optimal and Non-Pareto designs 




Table 4.4 shows the input as well as output of the runs that were carried out till including 
the last iteration. 

From the fig. 4.9, we can easily find out the Pareto solutions. The Pareto optimal points 
are marked and are shown in the fig. 4.10. And the hst of Pareto optimal points, stating 
their inputs as well as outputs is shown in Table 4.5. 

From Table 4.5, fig. 4.10 and fig. 4.1 1, it can be inferred that the Pareto optimal solutions 
are mainly consisting of the angles of attack around 30° and 37.5°. And the optimal points 
are different combinations of these angles with length of winglet starting fi-om 2.5 to 


3.75. 


1 Design Inputs 

Design Outputs 

Design 

No. 

Angle of Attack 

Length of Winglet 

Nusselt Number 

Pressure Drop 

1 

15 

2.5 

13.37 

0.0103 

2 

22.5 

2.5 

14.66 

0.015 

3* 

30 

2.5 

15.35 

0.016 

4' 

37.5 

2.5 

15.68 

0.017 

5' 

40 

2.5 

15.85 

0.02 

6* 

45 

2.5 

15.98 

0.023 

T 

50 

2.5 

16.23 

0.029 

8 

60 

2.5 

16.31 

0.038 

9 

15 

3.75 

14.31 

0.013 

10 

22.5 

3.75 

15.45 

0.019 

11' 

30 

3.75 

16.76 

0.035 



37.5 

3.75 

17.34 

0.038 

13' 

40 

3.75 

17.65 

0.04 

14' 

45 ' 

3.75 

18.2 

0.048 

15 

50 

3.75 

17.98 

0.054 

16 

30 

1.25 

12.34 

0.011 

17 

37.5 

1.25 

12.49 ^ 

0.0125 

18 

50 

1.25 

12.84 

0.015 

19' 

37.5 

3.125 

16.35 

0.03 

20* 

40 

3 125 

16.61 

0.034 

21 

45 

3.125 

16.75 

0.038 

22 

50 

3.125 

17.27 

0.045 

23' 

37.5 

3.5 

17.02 

0. 035 

24 

40 

3.5 

17.1 

0.038 

25 

45 

3.5 

17.42 

0 045 

26 

50 

3.5 

17.78 

0.053 

jiLKJ 

27 — 

37.5 

3.25 

16.41 

0.032 

28 

30 

2.25 

14.9 

0.0163 


Table 4.4; Input and Output of all the Designs 
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In Table 4.4 , the Pareto solutions are marked separately. 


Design Inputs 

Design Outputs 

Design 

No. 

Angle of Attack 

Length of Winglet 

Nusselt Number 

Pressure Drop 

1 

30 

2.5 

15.35 

0.016 

2 

37.5 

2.5 

15.68 

0.017 

3 

40 

2.5 

15.85 

0.02 

4 

45 

2.5 

15.98 

0.023 

5 

50 

2.5 

16.23 

0.029 

6 

30 

3.75 

16.76 

0.035 

7 

37.5 

3.75 

17.34 

0.038 

8 

40 

3.75 

17.65 1 

0.04 

9 

45 

3.75 

18.2 

0.048 

10 

37.5 

3.125 

16.35 

0.03 

11 

40 

3.125 

16.61 

0.034 

12 

37.5 

3.5 

17.02 

0.035 

13 

37.5 

3.25 

16.41 

0.032 


Table 4.5: Input and Output of Pareto Optimal Winglet Designs 
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5. Conclusions 


5.1 Technical Summary 

5.1.1 Results of the Present Work 

This study has appraised the utility of EVOP methods for deriving multiobjective 
solutions to engineering design problems, where GA or other methods cannot work due 
to practical limitations. This work illustrates that EVOP can be effectively used in such 
situations to find Pareto Optimal multi-objective engineering designs. 

Table 5.1 displays example Pareto Optiinal designs that were reached when optimization of channel 
geometry was attempted using the EVOP technique. These solutions shown possess Nusselt 
number between 15 and 18 roughly, while their pressure drop varies from 0.016 to 0.05. 


Design Inputs 

Design Outputs | 

Design 

No. 

Angle of Attack 

Length of Winglet 

Nusselt Number 

Pressure Drop 

1 

30 

2.5 

15.35 

0.016 

2 

37.5 

2.5 

15.68 

0.017 

3 

40 

2.5 

15.85 

0.02 

4 

45 

2.5 

15.98 

0.023 

5 

50 

2.5 

16.23 

0.029 

6 

30 

3.75 

16.76 

0.035 

7 

37.5 

3.75 

17.34 

0.038 

8 

40 ^ 

3.75 

17.65 

0.04 

9 

45 

3.75 

18.2 

0.048 

10 

37.5 

3.125 

16.35 

0.03 

11 

40 

3.125 

16.61 


12 

37.5 

3.5 

17.02 


13 

37.5 

3.25 

16.41 

0.032 


Table 5.1: Final Set of Pareto Optimal Winglet Designs 


5.1.2 Technical Conclusions 

The technical inferences specific to the heat exchanger design studied, of the present 
work may be summarized as follows. 
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1 . The angle of attack up to 22.5° does not provide us with considerable increment in 
Nusselt number as compared to larger angles of attack, though the pressure drop up to 
22.5° angle of attack is much lower comparatively. 

2. If we go for angles of attack larger than 40° then these designs show a large increment 
in pressure drop though there is increment in Nusselt number also. 

3. If we go on increasing angle of attack, then at a certain point the Nusselt number 
reaches its maximum, and if the angle of attack is forther increased beyond this point 
then it drops, and it is accompanied with larger pressure drop over previous angles of 
attack. 

4. In the present study, it is observed that the angle of attack between 30° and 40° give 
better results with Nusselt number in between 15 to 18, with a pressure drop between 
0.016 to 0.05 for the length of winglet between 2.5 to 3.75 

5. If the length of winglet is decreased below 2.5 then although the pressure drop is 
much lower, the Nusselt number does not increase much. So, at a lower length than 
2.5, there is a larger compensation for Nusselt number though pressure drop is 
improved somewhat. 

6. If the length of winglet is increased beyond 3.5, the increment in Nusselt number is 
always accompanied with larger and larger pressure drops. 

7. EVOP is a workable optimization technique for channel geometry decision making, 
as it can be used, in the present case, for calculating Nusselt number and Pressure 
drop in advance, and designmg of the channel and obstacle body will be much more 
convenient. It is very helpful in the sense that there is scarcity of knowledge as to 
what effects do the obstacle geometry has on the heat transfer phenomena. 


5.2 Overall Conclusions of the Present Work 

5.2.1 Optimization Methods and EVOP 

Evolutionary Operation has sometimes been referred to as an optimization technique and 
is compared with other optimization techniques. Such comparisons must be considered 
with some caution, for the following reasons. 


64 



Evolutionary Operation is an experimental technique for seeking the preferable. By 
contrast, problems of mathematical optimization are usually concerned with determining 
the maximum or minimum of some function z = / within some region R 

(xi, ,X]J of the space of the x’s, where usually it is assumed that: 

1 . The variables xi, ,Xk are known. 

2. The region R (xi, ,xi^ is known. 

3. The nature of the function z =/ ixi,....,x0 is known. 

4. The function z =/ {xi,....,xt^ can be computed without error for any chosen set of x’s. 
The situation which is met when we attempt to improve any process using EVOP differs 
from the mathematical optimization setup in a number of respects: 

1 . We do not know the variables xj,... .,Xk which should be included in the function/ 

2. We do not know with exactitude the region R (xj, ,X]^ in which we should seek to 

maximize z. 

3 . We do not know the functional form / 

4. The observations are usually subject to moderate or large errors. 

5.2.2 Consequences of using EVOP 

We may know, at any given stage of an EVOP investigation, that we have found 
conditions which are preferable to those previously considered. However we can rarely 
know that we have attained best conditions. Some of the reasons are as follows: 

1. Even if a process existed where improvement of a single response such as yield was 
the only criterion, we would seldom know how far or short of this theoretical point 
was the practical maximum reached by EVOP. 

2. In practice, the objective function will almost always need to take into account a 
number of responses such as yield, cost, and so on. It would be even less possible to 
be sure that the optimum objective had been attained when improvement was so 
many faceted. The Pareto concept may address static situation reasonably. However 
the importance of the various characteristics may change from time to time depending 
on the state of the market, the nature of competitive products, and so on. It is 
therefore extremely unlikely that we could ever be sure that the ultimate for a 
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multidimensional and shiftmg objective function of this kind had ever been achieved 
by EVOP. 


Most importantly, in the typical EVOP situation we cannot say at any point that no new 
possibilities for improvement will become available. Thus, optimum to EVOP is like the 
pot of gold at the end of a rainbow. 


5.3 Scope for Future Work 

1. The present study has used MAC algorithm and it has been limited to numerical 
investigation only. By execution of the physical heat transfer process, the outcome of 
this present work may be validated. This will require an experimental setup and the 
actual experimentation with different geometries. 

2. Some alternate, faster method than MAC algorithm, perhaps not as accurate, may be 
evaluated to broadly search the design space by EVOP before fine tuning by MAC 
algorithm is done. 

3 . The work is based on the EVOP technique. There are some methods like ROVOP 
(Rotating square Evolutionary Operation) and REVOP (Random Evolutionary 
Operations) and Simplex EVOP, which can also be used for the optimization of the 
same problem ( Black and Draper (1969)). A comparative study may show which 
method is more efficient in optimizing the heat exchanger design geometry. This will 
reduce the efforts and time in the future for bi-criteria heat exchanger geometry 
design. 

4 . The present study has taken into account only two design variables viz. angle of 
attack and the length of winglet obstacle. In order to help the decision making of heat 
exchanger design, aspect ratio, i.e. the ratio of width of channel to the height of 
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chaimel, also can be taken as the third input variable. This will uncover the effect of 
aspect ratio on heat enhancement and pressure penalty. This may be more helpful to 
us and lead to much better heat exchanger geometries. 


67 



References 

1. Adulbhan, P. and M. T. Tabucanon (1980). “Multicritenon Optimization in Industrial 
Systems Chapter 9, Decision models for industrial Systems Engineer and Managers, 
AIT, Bangkok, Thailand. 

2. Bagchi, T. P. (1999). Multiobjective Scheduling by Genetic Algorithms, Klluwer 

3. Box, G. E. P. and Draper, N. R. (1969), '^Evolutionary Operation: A Statistical 
Method for Process Improvement", John Wiley and Sons, Inc. New York. 

4. Cochran, J. L. And M. Zeleny, eds, (1973) "Multiple Criteria Decision Making”, 
University of South Carolina Press. 

5. Biswas, G., K. Torii, D. Fujii, K. Nishino "Numerical and experimental 
determination of flow structure and heat transfer effects of longitudinal vortices in a 
channel flow” , Int. J. Heat Transfer, Vol. 39, No. 16, pp. 3441-3451, 1996. 

6. Biswas, G., (1995), Solution of Navier-Stokes Equations for Incompressible Flows 
Using MAC and SIMPLE Algorithms, in Computational Fluid Flow and Heat 
Transfer, ed. K. Muralidhar and T. Sundararajan, Narosa Publishing House, India. 

7. Deb, K. And Srimvas, N(1995), "Multiobjective Optimization using Nondimensional 
Sorting Genetic Algorithm " MIT Press, Evolutionary Computations, 2(3), 221-248. 


8. Holman, J. P., "Heat Transfer”. McGraw-Hill Book Company, Singapore, 1989. 

9. Deb, Prashanta, (April, 1994), "Three Dimensional Computation of Laminar and 
Turbulent Flows with Heat Transfer in a channel with Built-in Longitudinal Vortex 
Generators ”, PhD. Thesis, ME Department, HT Kanpur. 


68 



10. Rao, P. Narayan (1999), “Multiobjective Optimization of Injection Molding Process 

Parameters: A Study using Genetic Algorithms”, M Tech. Thesis, IME, IIT Kanpur. 

■» * 

ll.Srinivas, T. D. (1998) “Multiple Criteria Optimization: An Evolution of Two 
Evolutionary Meta-Heuristic Methods ”, M Tech. Thesis, IME, IIT Kanpur. 

12. Tabucanon, M. T. (1989) “Multicriteria Decision Making in Industry”, Elsevier 
Science Publishers, New York. 


69 



